
#   APPENDIX
#   Figure 4    Differences in Marginal Effects of Price by Marginalization Threshold
#   Figure 5    Differences in Marginal Effects of Price by Poverty Threshold
#   Figure 6    Middle Class Messaging
#   Table 2     Flattening Participation: Pairwise Tests
#   Table 3      Sensitivity Analysis: Price Effects

# NEW TABLES HERE

#   Figure 7    Flattening Representation
#   Table 5     Heterogenous Effects: Voters
#   Table 6     Heterogenous Effects: Constituency and MP characteristics
#   Figure 9   ICT Use and Ownership in Africa
#   Figure 10   ICT Use and Ownership: Selected African Countries
#   Figure 11   ICT Use and Ownership Across Regions
#   Table 8     Global Inequality




#######################################################################################################
#                           APPENDIX
#######################################################################################################

#######################################################################################################
#   Figure 4 & 5    Difference in Marginal Effect of Price across all thresholds
#######################################################################################################

## exclude obs for which we are missing wealth info, there is no missingness in marginalization var
D1 <- DS[which(is.na(DS$POOR) == F), ] 
PM1 <- PM[which(is.na(DS$POOR) == F), ] 


vecWEALTH <- seq(0, .99, .01) #vector of population shares .01, .02
populate <- rep(NA, 100)   
populate <- quantile(D1$POOR, c(vecWEALTH), na.rm = TRUE)  #looping over all possible population splits
WEALTHdums <- function(x){
  (D1$POOR <x)
}
WEALTHMAT <- sapply(populate, WEALTHdums) #create matrix where columns are WEALTHY dummies at each level
WEALTHMAT <- WEALTHMAT*1
  
  
  
vecmarg <- seq(0, .99, .01) #vector of populaion shares .01, .02
populate <- rep(NA, 100)   
populate <- quantile(DS$MARG, c(vecmarg), na.rm = TRUE)  #looping over all possible population splits
margdums <- function(x){
  (DS$MARG <x)
}
MARGMAT <- sapply(populate, margdums) #create matrix where columns are marginal dummies at each level
MARGMAT <- MARGMAT*1

X = matrix(NA, 100, 4)


for (d in 2:100){
  wealth_res <- ri.inter.nocontrol(Y=D1$SMS, X=WEALTHMAT[,d], T=D1$PRICE, BLOCK=D1$BLOCK, PM=PM1)  #use D1 and PM1 to take into account missingness
  X[d,1:2] <- wealth_res[,1]
  marg_res <- ri.inter.nocontrol(Y=DS$SMS, X=MARGMAT[,d], T=DS$PRICE, BLOCK=DS$BLOCK, PM=PM)
  X[d,3:4] <- marg_res[,1]
}

save(X, file = "Results_Matrix")


#X[,2] <- substring(X[,2], 2, nchar(X[,2])-1)
#X[,4] <- substring(X[,4], 2, nchar(X[,4])-1)

pdf(file="../L/figures/MarginalEffectinequality.pdf", height=6, width=6)
plot(0,0, xlim= c(0,1), ylim=c(-.15,1), xlab="Threshold for classification as `poor'", ylab="Difference in marginal effect of subsidy", type="n")
points(vecWEALTH,X[,1], type="l", col="black", width=3)
points(vecWEALTH,X[,2], type="l", col="grey")
text(.82, 0.24, ("p-value"), col="grey") 
abline(0,0)
abline(h=.05, col= "red")
abline(h=.1, col= "red")
axis(4, at=c(.05, .1),
     col.axis="red", las=2, cex.axis=0.5, tck=-.01)
dev.off()



pdf(file="../L/figures/MarginalEffectMarginal.pdf", height=6, width=6)
plot(0,0, xlim= c(0,1), ylim=c(-.15,1), xlab="Threshold for classification as `marginalized'", ylab="Difference in marginal effect of subsidy", type="n")
points(vecmarg,X[,3], type="l", col="black")
points(vecmarg,X[,4], type="l", col="grey")
text(.82, 0.20, ("p-value"), col="grey") 
abline(0,0)
abline(h=.05, col= "red")
abline(h=.1, col= "red")
axis(4, at=c(.05, .1),
     col.axis="red", las=2, cex.axis=0.5, tck=-.01)
dev.off()

####################################################################################### 
# FIGURE 6: Middle Class Messaging
#######################################################################################


# logit regression of SMS on  wealth rank and wealth ranks squared
M     = glm(DS$SMS~DS$WEALTH_RANK + I(DS$WEALTH_RANK^2), data=DS, family = "binomial")
summary(M)
pred  <- predict(M,se.fit = TRUE)
const <- qt(0.975, df = M$df.residual)

#  bands on scale of link function
bands <- data.frame(X=DS$WEALTH_RANK[!is.na(DS$WEALTH_RANK)])
bands$upper <- pred$fit + (const * pred$se.fit)
bands$pred  <- pred$fit 
bands$lower <- pred$fit - (const * pred$se.fit)

# transform to response scale
bands <- data.frame(lapply(bands, binomial(link = "logit")$linkinv))

# plot
pdf(file="..L/figures/MC.pdf", width=8, height=6) 
plot(DS$WEALTH_RANK[!is.na(DS$WEALTH_RANK)], bands$pred, ylim = c(0, .07), xlab = "Wealth (Percentile)", ylab = "Probability of sending an SMS message", cex=.2)
points(DS$WEALTH_RANK[!is.na(DS$WEALTH_RANK)], bands$upper, col = "red", cex=.2)
points(DS$WEALTH_RANK[!is.na(DS$WEALTH_RANK)], bands$lower, col = "red", cex=.2)
dev.off()



#######################################################################################################
#   Table 1    Flattening Participation: Pairwise Tests
#######################################################################################################


# Table 1 appendix
xtable(h1tabf(D$V6_SMS==1, Y=D$POORBI, pctile = 1- mean(D$SMS[D$V6_SMS==1]), YLAB = "Share of poor respondents"))
xtable(h1tabf(D$V6_SMS==1, Y=D$WOMEN, pctile = 1- mean(D$SMS[D$V6_SMS==1]), YLAB = "Share of women respondents"))
xtable(h1tabf(D$V6_SMS==1, Y=D$nonCOE, pctile = 1- mean(D$SMS[D$V6_SMS==1]), YLAB = "Share of  non-coethnic respondents"))
xtable(h1tabf(D$V6_SMS==1, Y=D$nonCOG, pctile = 1- mean(D$SMS[D$V6_SMS==1]), YLAB = "Share of non-cogender respondents"))
xtable(h1tabf(D$V6_SMS==1, Y=D$DISTANT, pctile = 1- mean(D$SMS[D$V6_SMS==1]), YLAB = "Share of distant respondents"))


####################################################################################### 
# Table 2     Sensitivity analysis: Price Effects on Message Content
####################################################################################### 
# recodefine pub/priv var and recreate table 6
DS$SMS_Priv <- rep(0, 3790)
DS$SMS_Priv[DS$PubPriv==0 | DS$PubPriv==1 | DS$PubPriv==2 |DS$PubPriv==4] <- 1

DS$SMS_Pub <- rep(0, 3790)
DS$SMS_Pub[DS$PubPriv==3 | DS$PubPriv==5] <- 1


#top part of the table
MainT=cbind(
  sapply(0:2, function(i) mean(DS$SMS[DS$PRICE==i])),
  sapply(0:2, function(i) mean(DS$SMS_Pub[DS$PRICE==i])),
  sapply(0:2, function(i) mean(DS$SMS_Priv[DS$PRICE==i])),
  sapply(0:2, function(i) sum(DS$PRICE==i, na.rm=T)))

rownames(MainT) <- c("Full", "Subsidy", "Free")
colnames(MainT) <- c("Any SMS", "Public SMS", "Private SMS", "Obs")
xtable(MainT, digits=3)

# BOTTOM PART OF THE TABLE

# DEFINE REDUCED Y
Y_REDU2=cbind(DS$SMS, DS$SMS_B, DS$SMS_Priv)
OUT = matrix(NA, 12,3)
for(y in 1:3){
  OUT[1:3,y]     <-ri_lm(Y_REDU2[,y], (DS$PRICE==1), DS$BLOCK, PM, C = (DS$PRICE!=2), addn=TRUE, round=3, onesided=3)  # T_FS  
  OUT[4:6,y]     <-ri_lm(Y_REDU2[,y], (DS$PRICE==2), DS$BLOCK, PM, C = (DS$PRICE!=0), addn=TRUE, round=3, onesided=3)  # T_SO
  OUT[7:9,y]     <-ri_lm(Y_REDU2[,y], (DS$PRICE==2), DS$BLOCK, PM, C = (DS$PRICE!=1), addn=TRUE, round=3, onesided=3)  # T_FO
}

for(y in 1:3){
  OUT[10:12,y]   <-ri_lm(Y_REDU2[,y],  DS$PRICE,  DS$BLOCK, PM, addn=TRUE, round=3, onesided=3)   # Linear   
}

# H4 test using SUR, Appendix version
msurH4A 	<- systemfit(list(MV = SMS_Priv~PRICE + as.factor(BLOCK), MB = SMS_Pub~PRICE + as.factor(BLOCK)), data=DS)
CM_H4A	=  coef(summary(msurH4A))
LH_H4A 	= linearHypothesis(msurH4A, "MB_PRICE  - MV_PRICE", test = "Chisq", na.rm=TRUE)
pH4_plusA = round(LH_H4A[[4]][2]/2,3)
pH4_plusA

OUTA= rbind(cbind(OUT, c(rep("",9),pH4_plusA)), c("H2 test","","",""))
C1A = c("Subsidy vs. Full Price", "", "","Free vs. Subsidy","", "","Free vs. Full Price","", "", "Linear Trend", "","", "")
C2A = c(rep(c("ATE", paste("(p)"), "(N)"),3), c("Trend", "(p)", "(N)"), "")
OUT2A = cbind(C1A,C2A, OUTA)
colnames(OUT2A)<-c("Treatment", "Effect", "Any", "Public","Private", "H4 test")
print(xtable(OUT2A, digits=3), include.rownames=FALSE)




#################################################################################################################
# # Table 3  
# Reproduce Table 6 from main analysis in Appendix (Price Induced Flattening 1: Test of H 3.1)  using RI
#################################################################################################################


OUT9 <- matrix(NA, 6, 4)

#Rich
OUT9[2:3,1] <- as.matrix(ri_sg(Y=DS$SMS, X=DS$POORBI, T=DS$PRICE, BLOCK=DS$BLOCK, control=data.frame(DS$ACC), PERM=PM, effecttype=2))
#Poor
OUT9[2:3,2] <- as.matrix(ri_sg(Y=DS$SMS, X=DS$POORBI, T=DS$PRICE, BLOCK=DS$BLOCK, control=data.frame(DS$ACC), PERM=PM, effecttype=1))

#Low Access
OUT9[5:6,1] <- as.matrix(ri_sg(Y=DS$SMS, X=DS$ACC, T=DS$PRICE, BLOCK=DS$BLOCK, control=data.frame(DS$POORBI), PERM=PM, effecttype=2)) 
#High Access
OUT9[5:6,2] <-  as.matrix(ri_sg(Y=DS$SMS, X=DS$ACC, T=DS$PRICE, BLOCK=DS$BLOCK, control=data.frame(DS$POORBI), PERM=PM, effecttype=1))


#########Marginal effect of subsidy by POOR controlling for access

OUT9[2:3,4]  <- ri.inter(Y=DS$SMS, X=DS$POORBI, T=DS$PRICE, control=data.frame(DS$ACC), BLOCK=DS$BLOCK, PM=PM)

#########Marginal effect of subsidy by ACCESS controlling for poor

OUT9[5:6,4]  <-  ri.inter(Y=DS$SMS, X=DS$ACC, T=DS$PRICE, control=data.frame(DS$POORBI), BLOCK=DS$BLOCK, PM=PM)

rows = c("", "Marginal Effect of Subsidy by Poverty","","", "Marginal Effect of Subsidy by Political Access", "")
OUT9 <- cbind(rows, OUT9)
OUT9[1, 1:5] <- c("", "Rich", "Poor", "Hypothesized Difference", "Above-Below")
OUT9[4, 1:5] <- c("", "Low Access", "High Access", "Hypothesized Difference", "Above-Below")

xtable(OUT9)


#######################################################################################################
# Table 4
# Reproduce Table 7 from main analysis in Appendix (Price Induced Flattening 2: Test of H 3.2)  using RI
#######################################################################################################
OUT7 <- matrix(NA, 3 ,2)

#Marginalized
OUT7[1,1:2] <- as.matrix(ri_sg(Y=DS$SMS, X=DS$MARGINAL, T=DS$PRICE, BLOCK=DS$BLOCK, control=NULL, PERM=PM, effecttype=1))

#Not Marginalized
OUT7[2,1:2]  <- as.matrix(ri_sg(Y=DS$SMS, X=DS$MARGINAL, T=DS$PRICE, BLOCK=DS$BLOCK, control=NULL, PERM=PM, effecttype=2))

# Difference Marginalized v Non-Marginalized, no controls
OUT7[3,1:2] <- ri.inter.nocontrol(Y=DS$SMS, X=DS$MARGINAL, T=DS$PRICE, BLOCK=DS$BLOCK, PM=PM)


OUT7 <- round(OUT7, 4)

OUT7 <- cbind(c("Marginal Effect of Price on Marginalized", "Marginal Effect of Price on Non-marginalized", "Difference"), OUT7)
xtable(OUT7)

SEND = DS$SMS==1
h(DS$MARGINAL[SEND], DS$PRICE[SEND])
bss=bsh(DS$MARGINAL[SEND],PM[SEND,])
hist(bss)

# changed to reflect treatment as subsidy (0 is full price, 1 is subsidy, 2 is complete subsidy)
h5.3tab = function(Y, X, name, PERM=PM[SEND,]) (c(
  paste("Share of ",  name, " among high price senders=", round(mean(Y[X==0]),3), sep=""), 
  paste("n_0  =", round(length(Y[X==0]),3)),
  paste("Share of ",  name, " among mid price senders= ",round(mean(Y[X==1]),3), sep=""),
  paste("n_S  =", round(length(Y[X==1]),3)),
  paste("Share of ",  name, " among low price senders= ",round(mean(Y[X==2]),3), sep=""),
  paste("n_F  =", round(length(Y[X==2]),3)),
  paste("Trend from high to low=",  round(h(Y, X),3)), 
  paste("(p (neg))=",mean(h(Y, X)>=bsh(Y, PERM))),
  paste("(p (two sided))=",mean(abs(h(Y, X))>=abs(bsh(Y, PERM)))),
  paste("(p (pos))=",mean(h(Y, X)<=bsh(Y, PERM)))
))
h5.3tab(DS$MARGINAL[SEND],DS$PRICE[SEND], "marginal respondents")
#h5.3tab(DS$POORBI[SEND],DS$PRICE[SEND], "poor respondents")
#h5.3tab(DS$WOMEN[SEND],DS$PRICE[SEND], "women")




#######################################################################################################
# Figure 7  Flattening Representation 
#######################################################################################################
# Function to construct test (change breaks to 0:13) 
Dev.Dist.Out  = function(D, sims, rond=3, 
                         xyxlab="SMS user deviation", xyylab="Engaged deviation",  
						 breaks=0:13, cex=.2, xlim=NULL, ylim=NULL, breaks2=50){
  s = as.vector(table(D$SMS,D$ENG)) # (00), (SMS=1, ENG=0), (SMS=0, ENG=1), (11)
  s = s/(sum(s))           # Shares in each group
  
  # Function for prediction from Multinomial logit
  mlpred	= function(SMS, ENG, B){(exp(B)%*%c(1,SMS,ENG))/sum(exp(B)%*%c(1,SMS,ENG))}
  # Prediction of shares for a given type	
  mlpred_b	= function(s,B){s[1]*mlpred(0,0,B)+s[2]*mlpred(1,0, B)+s[3]*mlpred(0,1,B)+s[4]*mlpred(1,1,B)}  
  # Prediction given the different types that exist in the population; weighted by share
  
  DEVIATION = function(SMS,ENG, B){
		sum((mlpred(SMS,ENG, B) - mlpred_b(s,B))^2)/2
		}
  # NRS statistic for difference betwee a particular population(SMS,ENG combination) and the full population; given B		
		
  # Run DEVIATION for matrix of simulated betas
  Dev.Dist = function(SMS, ENG, BS,  sims){ 
		sapply(1:sims, function(i) {mu = t(matrix(BS[i,],3)); DEVIATION(SMS,ENG, mu)})
		}
  
  # Multinomial Model
  #############################################################################
  M <- multinom(priority ~ SMS + ENG, data = D)
  # Simulated betas
  BS = mvrnorm(sims, mu = as.vector(t(coef(M))), Sigma = vcov(M))
  BS = cbind(matrix(0, sims, 3), BS)
  # Quantity of interest
  SMSDist = Dev.Dist(1,0, BS, sims=sims)
  ENGDist = Dev.Dist(0,1, BS, sims=sims)
  DiffDist = ENGDist - SMSDist  # Extent to which ENG (engaged) dist is large

  par(mfrow=c(2,2))
	  opar=par(ps=10)
	  
	  hist(SMSDist,   main = paste("SMS (only) group deviation:", 
								   round(mean(SMSDist),rond)), xlab="Deviation", xlim=c(0, max(c(SMSDist,ENGDist))))
	  hist(ENGDist,   main = paste("Engaged (only) group deviation:", 
								   round(mean(ENGDist),rond)), xlab="Deviation", xlim=c(0, max(c(ENGDist))))
	  plot(SMSDist[SMSDist>ENGDist], ENGDist[SMSDist>ENGDist], ylab = xyylab, xlab = xyxlab, main="Simulated deviations",
		   xlim=xlim,
		   ylim=ylim,
		   cex=cex, col="blue"); 
	  points(SMSDist[SMSDist<=ENGDist], ENGDist[SMSDist<=ENGDist],col="red", cex=cex)
	  abline(a=0, b=1, col="black"); 
	  points(mean(SMSDist), mean(ENGDist),col="black", pch=16); 
	  hist(DiffDist, main = paste("Distribution of the difference in deviations  \n Mean=", 
								  round(mean(DiffDist),rond),  ", sd=", 
								  round(sd(DiffDist),rond),   ", % negative:",  round(mean(DiffDist<0),rond), sep=""), 
		   xlab="Difference in Deviation", breaks=breaks2)
	  abline(v=0, col="red")
	}


# Graph
pdf(file="../L/figures/NonRepRealData.pdf", width=8, height=9) 
	Dev.Dist.Out(D, 5000, cex=.2, xlim = c(0,.09), ylim=c(0, .09))
dev.off()



#######################################################################################################
#   Table 6   Heterogeneous Effects: Voters' Characteristics (expanded set of results)
#######################################################################################################
# Row and Column Conditions for Table
rC = cbind((DS$ACC==0), (DS$ACC==1), TRUE)
cC = cbind((DS$POORBI==0), (DS$POORBI==1), TRUE)
DV = cbind(DS$SMS_V, DS$SMS_B, DS$SMS)

# Output
X = matrix(NA, 24,4)

for(d in 1:3)for(r in 1:3)for(c in 1:3){X[((r-1)*6+(d-1)*2+1):((r-1)*6+(d-1)*2+2), c]<- {
  ri_lm(DV[,d], DS$PRICE, DS$BLOCK, PM=PM,  C=(rC[,r]&cC[,c]), onesided=3)}} 


# 1 private messaging outcome, predict interaction on POORBI given low access
C = (DS$ACC == 0)
X[1:2,4] <- ri.inter.nocontrol(Y=DS$SMS_V[C], X=DS$POORBI[C], T=DS$PRICE[C], BLOCK=DS$BLOCK[C], PM=PM[C,])

# 2 public messaging outcome, predict interaction on POORBI given low access
X[3:4,4] <- ri.inter.nocontrol(Y=DS$SMS_B[C], X=DS$POORBI[C], T=DS$PRICE[C], BLOCK=DS$BLOCK[C], PM=PM[C,])

# 3 any messaging outcome, predict interaction on POORBI given low access
X[5:6,4] <- ri.inter.nocontrol(Y=DS$SMS[C], X=DS$POORBI[C], T=DS$PRICE[C], BLOCK=DS$BLOCK[C], PM=PM[C,])


# 4 private messaging outcome, predict interaction on POORBI given high access
B = (DS$ACC == 1 )
X[7:8,4]<- ri.inter.nocontrol(Y=DS$SMS_V[B], X=DS$POORBI[B], T=DS$PRICE[B], BLOCK=DS$BLOCK[B], PM=PM[B,])
# 5 public messaging outcome, predict interaction on POORBI given high access
X[9:10,4]<- ri.inter.nocontrol(Y=DS$SMS_B[B], X=DS$POORBI[B], T=DS$PRICE[B], BLOCK=DS$BLOCK[B], PM=PM[B,])
# 6 any messaging outcome, predict interaction on POORBI given high access
X[11:12,4]<- ri.inter.nocontrol(Y=DS$SMS[B], X=DS$POORBI[B], T=DS$PRICE[B], BLOCK=DS$BLOCK[B], PM=PM[B,])
# 
#  NO SUBSETTING
# 7 private messaging outcome, predict interaction on POORBI given any access
X[13:14,4]<- ri.inter.nocontrol(Y=DS$SMS_V, X=DS$POORBI, T=DS$PRICE, BLOCK=DS$BLOCK, PM=PM)
# 8 public messaging outcome, predict interaction on POORBI given any access
X[15:16,4]<- ri.inter.nocontrol(Y=DS$SMS_B, X=DS$POORBI, T=DS$PRICE, BLOCK=DS$BLOCK, PM=PM)
# 9 any messaging outcome, predict interaction on POORBI given any access
X[17:18,4]<- ri.inter.nocontrol(Y=DS$SMS, X=DS$POORBI, T=DS$PRICE, BLOCK=DS$BLOCK, PM=PM)


# ###Lower rows column 1
# 1 private messaging outcome, interaction on access given rich
G = (DS$POORBI == 0)
X[19:20,1]<-  ri.inter.nocontrol(Y=DS$SMS_V[G], X=DS$ACC[G], T=DS$PRICE[G], BLOCK=DS$BLOCK[G], PM=PM[G,])
# 2 public messaging outcome, interaction on access given rich
X[21:22,1]<- ri.inter.nocontrol(Y=DS$SMS_B[G], X=DS$ACC[G], T=DS$PRICE[G], BLOCK=DS$BLOCK[G], PM=PM[G,])
# 3 any  messaging outcome, interaction on access given rich
X[23:24,1]<-  ri.inter.nocontrol(Y=DS$SMS[G], X=DS$ACC[G], T=DS$PRICE[G], BLOCK=DS$BLOCK[G], PM=PM[G,])

###Lower rows column 1
# 4 private messaging outcome, interaction on access given poor
E = (DS$POORBI == 1)
X[19:20,2]<- ri.inter.nocontrol(Y=DS$SMS_V[E], X=DS$ACC[E], T=DS$PRICE[E], BLOCK=DS$BLOCK[E], PM=PM[E,])
# 5 public messaging outcome, interaction on access given poor
X[21:22,2]<- ri.inter.nocontrol(Y=DS$SMS_B[E], X=DS$ACC[E], T=DS$PRICE[E], BLOCK=DS$BLOCK[E], PM=PM[E,])
# 6 any  messaging outcome, interaction on access given pooR
X[23:24,2]<- ri.inter.nocontrol(Y=DS$SMS[E], X=DS$ACC[E], T=DS$PRICE[E], BLOCK=DS$BLOCK[E], PM=PM[E,])

##Lower rows, column 3
# 7 private messaging outcome, interaction on access, given any wealth
X[19:20,3]<- ri.inter.nocontrol(Y=DS$SMS_V, X=DS$ACC, T=DS$PRICE, BLOCK=DS$BLOCK, PM=PM)
# 8 public messaging outcome, interaction on access, given any wealth
X[21:22,3]<- ri.inter.nocontrol(Y=DS$SMS_B, X=DS$ACC, T=DS$PRICE, BLOCK=DS$BLOCK, PM=PM)
# 9 any  messaging outcome, interaction on access, given any wealth
X[23:24,3]<- ri.inter.nocontrol(Y=DS$SMS, X=DS$ACC, T=DS$PRICE, BLOCK=DS$BLOCK, PM=PM)

class(X) <- "numeric"
X <- round(X, digits=4)

# Export Table
Xtest=cbind(X[,1:2], c(rep("",22), "H_3.1", "(+)"),  X[,3:ncol(X)])
Xtest= cbind(Xtest[,1:4], c(rep("",16), "H_3.2", "(+)", rep("", 6)), Xtest[,5:ncol(Xtest)])

labels=as.matrix(cbind(c("Low Access", rep("", 5), "High Access", rep("", 5),
          "Any Access", rep("", 5),"Difference", rep("", 5)),
          c(rep(c("Private","","Public","", "Any", ""),4))))

Xtest= cbind(labels, Xtest)
        
colnames(Xtest)<-c("", "", "Rich", "Poor", "", "All", "", "Difference")
print(xtable(Xtest),  include.rownames=FALSE)
# note: p values are ones sided (positive) for all effects





#######################################################################################################
#   Table 7  Constituency and MP Characterisitics
#######################################################################################################
options(digits=3)
HET= matrix(NA,12,3)

#uptake columns: for each COVAR (0,1) in the full price treatment, what % sent a message
HET[c(1),1] <- round(sum(DS$SMS[DS$NRM==0  & DS$PRICE==0])/(sum(DS$NRM==0 & DS$PRICE==0)), digits=3)
HET[c(1),2] <- round(sum(DS$SMS[DS$NRM==1  & DS$PRICE==0])/(sum(DS$NRM==1 & DS$PRICE==0)), digits=3)

HET[c(4),1]<- round(sum(DS$SMS[DS$COPART==0 & DS$PRICE==0], na.rm= TRUE)/sum((DS$COPART==0 &  DS$PRICE==0), na.rm= TRUE), digits=3)
HET[c(4),2]<- round(sum(DS$SMS[DS$COPART==1 & DS$PRICE==0], na.rm= TRUE)/sum((DS$COPART==1 &  DS$PRICE==0), na.rm= TRUE), digits=3)

HET[c(7),1] <-round(sum(DS$SMS[DS$COMPETITIVE== 0 &  DS$PRICE==0])/(sum(DS$COMPETITIVE== 0 &  DS$PRICE==0)), digits=3)
HET[c(7),2] <-round(sum(DS$SMS[DS$COMPETITIVE== 1 &  DS$PRICE==0])/(sum(DS$COMPETITIVE== 1 &  DS$PRICE==0)), digits=3)

HET[c(10),1] <- round(sum(DS$SMS[DS$YOUNG== 0 & DS$PRICE==0])/(sum(DS$YOUNG== 0 &  DS$PRICE==0)), digits=3)
HET[c(10),2] <- round(sum(DS$SMS[DS$YOUNG== 1 & DS$PRICE==0])/(sum(DS$YOUNG== 1 &  DS$PRICE==0)), digits=3)


## tests for heterogeneous effects-- marginal effect of price conditional on:
# (a) NRM (b)  Copartisanship (c) competitiveness (d) MP age
# Test that marginal effect of subsidy is postive for subgroups of interest


# effect on NRM
HET[2:3, 2] <- as.matrix(ri_sg(Y=DS$SMS, X=DS$NRM, 		T=DS$PRICE, BLOCK=DS$BLOCK, PERM=PM, effecttype=1))
#effect on non- NRM
HET[2:3, 1] <- as.matrix(ri_sg(Y=DS$SMS, X=DS$NRM, 		T=DS$PRICE, BLOCK=DS$BLOCK, PERM=PM, effecttype=2))

#effect on copart
HET[5:6, 2] <- as.matrix(ri_sg(Y=DS$SMS, X=DS$COPART, 	T=DS$PRICE, BLOCK=DS$BLOCK, PERM=PM, effecttype=1))
#effect on non-copart
HET[5:6, 1] <- as.matrix(ri_sg(Y=DS$SMS, X=DS$COPART, 	T=DS$PRICE, BLOCK=DS$BLOCK, PERM=PM, effecttype=2))

#effect on compete
HET[8:9, 2] <- as.matrix(ri_sg(Y=DS$SMS, X=DS$COMPETITIVE, 	T=DS$PRICE, BLOCK=DS$BLOCK, PERM=PM, effecttype=1))
#effect on non-compete
HET[8:9, 1] <- as.matrix(ri_sg(Y=DS$SMS, X=DS$COMPETITIVE, 	T=DS$PRICE, BLOCK=DS$BLOCK, PERM=PM, effecttype=2))

A = as.matrix(ri_sg(Y=DS$SMS, X=DS$COMPETITIVE, 			T=DS$PRICE, BLOCK=DS$BLOCK, PERM=PM, effecttype=1))


#effect on young
HET[11:12, 2] <- as.matrix(ri_sg(Y=DS$SMS, X=DS$YOUNG, 	T=DS$PRICE, BLOCK=DS$BLOCK, PERM=PM, effecttype=1))
#effect on non-young
HET[11:12, 1] <- as.matrix(ri_sg(Y=DS$SMS, X=DS$YOUNG, 	T=DS$PRICE, BLOCK=DS$BLOCK, PERM=PM, effecttype=2))

# RI for interaction effects  
HET[2:3, 3]  <- ri.inter.nocontrol(Y=DS$SMS,  X=DS$NRM, T=DS$PRICE, BLOCK=DS$BLOCK, PM=PM)
## replace NAs in copartisanship with 0s, if you did not give us data on partisanship we assume you are not aligned with your MP
DS$COPART[is.na(DS$COPART)] <- 0
HET[5:6, 3]  <- ri.inter.nocontrol(Y=DS$SMS,  X=DS$COPART, 		T=DS$PRICE, BLOCK=DS$BLOCK, PM=PM)
HET[8:9, 3]  <- ri.inter.nocontrol(Y=DS$SMS,  X=DS$COMPETITIVE, T=DS$PRICE, BLOCK=DS$BLOCK, PM=PM)
HET[11:12, 3] <- ri.inter.nocontrol(Y=DS$SMS, X=DS$YOUNG, 		T=DS$PRICE, BLOCK=DS$BLOCK, PM=PM)


C1 = c("NRM", "", "","Copartisan","", "","Competitive","", "", "Younger", "","")
C2 = c(rep(c("Uptake", "ATE", paste("(p)")),4))
OUT = cbind(C1,C2, HET)
colnames(OUT)<-c("","", "No", "Yes", "ATE Difference")
print(xtable(OUT, digits=3), include.rownames=FALSE)

# Checks
coef(summary(lm(SMS ~ COMPETITIVE*PRICE + as.factor(BLOCK), data  = DS)))[c(2:3,297),]
coef(summary(lm(SMS ~ YOUNG*PRICE + as.factor(BLOCK), data  = DS)))[c(2:3,297),]



############################################################################
#   Figure 9   Rates of ICT use/ownership of 10 years: all of SSA
############################################################################
ICT <- read.dta("data/ITU_dataset.dta")

ICT$countryid <- as.numeric(factor(with(ICT, paste(ICT$Country))))  # Unique country ID
ICT$regionid <- as.numeric(factor(with(ICT, paste(ICT$Region))))    # Unique region ID
all.years <- as.list(ICT$year)


Africa <- ICT[which(ICT$region==1), ]

#number of countries
colors <- length(unique(Africa$countryid))
countries <- unique(Africa$countryid)


# get the range for the x and y axis of each plot
xrange  <- range(Africa$year) 
yrange1 <- range(Africa$HHComputer, finite=TRUE) 
yrange2 <- range(Africa$HHInternetAccess, finite=TRUE)
yrange3 <- range(Africa$InternetUsers, finite=TRUE)
yrange4 <- range(Africa$MobileSubscriptions, finite=TRUE)
yrange <- rbind(yrange1, yrange2, yrange3, yrange4)

labels <- c("Percent HHs with Computer", "Percent HHs with Internet Access", "Percent Internet Users", "Mobile Subscribers per 100 people (Users)")
titles <- c("Computer ownership", "Internet Access", "Internet Use", "Mobile Subscriptions")
vars  <-  c("HHComputer", "HHInternetAccess", "InternetUsers", "MobileSubscriptions")

pdf("../L/figures/AfricaICT.pdf")
op <- par(
  oma=c(0,0,0,8), 
  mfrow=c(2,2)
)

colors <- rainbow(length(countries))

for (j in 1:4){
  plot(xrange, yrange[j, ], type="n", xlab="Years",
       ylab=paste(labels[j]) ) 
  
  title(paste(titles[j]) ) 
  
  for (i in 1:length(countries)){ 
    country <- subset(Africa, Africa$countryid==countries[i]) 
    lines(country$year, country[,vars[j]], type="o", lwd=.5, col=colors[i], pch=20)
  }
}

par(op)                      # Leave the last plot
op <- par(usr=c(0,1,0,1),    # Reset the coordinates
          xpd=NA)            # Allow plotting outside the plot region
legend('right', 
       unique(Africa$Country), cex=.5, pt.cex = .8, col=colors,
       pch=20, title="Countries",)

dev.off()


########################################################################################
### Figure 10 ICT Use and Ownership for Selected African Countries
########################################################################################


# IDs of selected African countries
IDs <- c(5, 30, 28, 42, 58, 66, 89, 108, 129, 148, 157, 182, 194, 195)

SmallAf <- subset(ICT, ICT$countryid==5 | ICT$countryid==30 | ICT$countryid==28 | ICT$countryid==42 | ICT$countryid==58 |
                    ICT$countryid==66 | ICT$countryid==89 | ICT$countryid==108 | ICT$countryid==129 | ICT$countryid==148 | ICT$countryid==157 |
                    ICT$countryid==182 | ICT$countryid==194 | ICT$countryid==195)


#number of countries
colors <- length(unique(SmallAf$countryid))
countries <- unique(SmallAf$countryid)


# get the range for the x and y axis of each plot
xrange  <- range(SmallAf$year) 
yrange1 <- range(SmallAf$HHComputer, finite=TRUE) 
yrange2 <- range(SmallAf$HHInternetAccess, finite=TRUE)
yrange3 <- range(SmallAf$InternetUsers, finite=TRUE)
yrange4 <- range(SmallAf$MobileSubscriptions, finite=TRUE)
yrange <- rbind(yrange1, yrange2, yrange3, yrange4)


pdf("../L/figures/SmallAfICT.pdf")
op <- par(
  oma=c(0,0,0,8), 
  mfrow=c(2,2)
)

colors <- rainbow(length(countries))

for (j in 1:4){
  plot(xrange, yrange[j, ], type="n", xlab="Years",
       ylab=paste(labels[j]) ) 
  
  title(paste(titles[j]) ) 
  
  for (i in 1:length(countries)){ 
    country <- subset(SmallAf, SmallAf$countryid==countries[i]) 
    lines(country$year, country[,vars[j]], type="o", lwd=.5, col=colors[i], pch=20)
  }
}

par(op)                      # Leave the last plot
op <- par(usr=c(0,1,0,1),    # Reset the coordinates
          xpd=NA)            # Allow plotting outside the plot region
legend('right', 
       unique(SmallAf$Country), cex=.5, pt.cex = .8, col=colors,
       pch=20, title="Countries",)

dev.off()


################################################################
# Figure 11 ICT Ownership across regions
################################################################

uganda <- subset(ICT, ICT$countryid==182) 


#number of Regions
regions <- unique(ICT$regionid)


# get the range for the x and y axis of each plot
xrange  <- range(ICT$year) 
yrange1 <- range(ICT$HHComputer, finite=TRUE) 
yrange2 <- range(ICT$HHInternetAccess, finite=TRUE)
yrange3 <- range(ICT$InternetUsers, finite=TRUE)
yrange4 <- range(ICT$MobileSubscriptions, finite=TRUE)
yrange <- rbind(yrange1, yrange2, yrange3, yrange4)


pdf("../L/figures/RegionICT.pdf")
op <- par(
  oma=c(0,0,0,8), 
  mfrow=c(2,2)
)

colors <- rainbow(length(regions))
ICTnew <- c(ICT$year, ICT$regionid, ICT$HHComputer, ICT$HHInternetAccess, ICT$InternetUsers, ICT$MobileSubscriptions)


for (j in 1:4){
  plot(xrange, yrange[j, ], type="n", xlab="Years",
       ylab=paste(labels[j]) ) 
  title(paste(titles[j]) ) 
  
  
  means <- as.data.frame(aggregate(ICT, by=list(ICT$year,ICT$regionid), FUN=mean, na.rm=TRUE))  
  
  for (i in 1:length(regions)){ 
    region <- subset(means, means$regionid==regions[i]) 
    lines(region$year, region[,vars[j]], type="o", lwd=.5, col=colors[i], pch=20)    
    lines(uganda$year, uganda[,vars[j]], type="o", lwd=.5, col="black", pch=20)
    
  }
}


par(op)                      # Leave the last plot
op <- par(usr=c(0,1,0,1),    # Reset the coordinates
          xpd=NA)            # Allow plotting outside the plot region
legend('right', 
       unique(ICT$Region), cex=.5, pt.cex = .8, col=colors,
       pch=20, title="Countries",)

dev.off()


#######################################################################################################
#  Table 8 Inequality
#######################################################################################################
GINI <- read.csv("data/WB_GINI.csv")
xtable(GINI)

#######################################################################################################
#   Summary Stats
#######################################################################################################
stargazer(D)

#######################################################################################################
# Figure XX  Private and Public Messaging by Group		
#######################################################################################################

se <- function(x) sqrt(var(x,na.rm=TRUE)/length(na.omit(x)))
extract.mean = function(X, val = 1){c(mean(DS$SMS_V[X==val], na.rm=TRUE), 	mean(DS$SMS_B[X==val], na.rm=TRUE), mean(DS$SMS[X==val], na.rm=TRUE))}
extract.se = function(X, val = 1){c(se(DS$SMS_V[X==val]), se(DS$SMS_B[X==val]), se(DS$SMS[X==val]))}
extract.mean(DS$nonCOE)
extract.se(DS$nonCOE)
rnames=c("Co-Ethnic", "Co-Gender", "Close to district capital", "Formal Education", "Wealthier", "Men", "Privileged (Less Marginalized)", "High Access", "Highly Engaged")

collect.mean = function(on=1){
  X= rbind( extract.mean(DS$nonCOE, 	  val=on),
            extract.mean(DS$nonCOG, 	  val=on),
            extract.mean(DS$DISTANT, 	  val=on),
            extract.mean(DS$EDUBI, 		  val=1-on),
            extract.mean(DS$POORBI, 	  val=on),
            extract.mean(DS$WOMEN, 		  val=on),
            extract.mean(DS$MARGINAL, 	  val=on),
            extract.mean(DS$ACC, 	  val=1-on),
            extract.mean(DS$ENGBI, val=1-on)
  )
  rownames(X) = rnames
  X	
}

collect.se = function(on=1){
  X= rbind( extract.se(DS$nonCOE, 	  val=on),
            extract.se(DS$nonCOG, 	  val=on),
            extract.se(DS$DISTANT, 	  val=on),
            extract.se(DS$EDUBI, 		  val=1-on),
            extract.se(DS$POORBI, 	  val=on),
            extract.se(DS$WOMEN, 		  val=on),
            extract.se(DS$MARGINAL, 	  val=on),
            extract.se(DS$ACC, 	  val=1-on),
            extract.se(DS$ENGBI, val=1-on)
  )
  rownames(X) = rnames
  X	
}
# First three columns are for marginalized, next three are for privileged
b 	= cbind(collect.mean(1), collect.mean(0))
se 	= cbind(collect.se(1), collect.se(0))

# b[1,1]<-0  # test  to make sure items are in correct order 
pdf(file="../L/figures/Senders.pdf", width=16, height=10) 
mat_plot2(b, se, title="", rnames = rnames,rnames_sub = c("No", "Yes"), cnames = c("Public", "Private", "Any"),   groups=2, 
          rnamestab=0, rset=.5, rset2=0.3, 
          norm=10, center=.04, vlinex=-.04, vcol="black", tick = .04, tick2=.02, ticktextsize=.8,
          bottom=1, uset=1, gap=.3)
dev.off()

